#################################
#Batista Pereira et al, JEPS
#Code for Appendix 6  - Randomization test
#Analysis for testing randomization in Study 1


#################################

library(tidyverse)
library(randomizr)
library(ggplot2)
library(xtable)
options(scipen=999) 

home <- "Replication/Randomization-Test" # set appropriate working directory here 

data_sim <- read_csv(paste0(home, "/quaest_radomsurvey_results.csv"))

#checking data
glimpse(data_sim)
table(data_sim$q2) #randomization question
table(data_sim$answer_button) #outcome 

#rename variable with treatment assignment
data_sim <- data_sim %>% 
  mutate(treat = ifelse(q2 == "Presidente Bolsonaro tem feito várias declarações contra a vacina. Qual a sua opinião sobre avacina: ela é segura ou insegura?", 1, 0))

#checking 
table(data_sim$treat)
table(is.na(data_sim$treat))

#simulating data
data_sim_R1 <- list() #
data_sim_R2 <- list()
set.seed(30047)
for (i in 1:1000){
  data_sim_R1[[i]] <- slice_sample(data_sim, n = 1000, replace = T) #creating a random sample of 1,000 respondents (1,000 times)
  data_sim_R2[[i]] <- simple_ra(N = nrow(data_sim_R1[[i]])) #random assignment for a sample of equal size to the empirical sample (1,000 times)
}

#Number of units assigned to treatment
func <- function(x) {
  sum(x$treat)
}

treats_quaest <- sapply(data_sim_R1, func)
treats_randomizr <- sapply(data_sim_R2, sum)

hist(treats_quaest)
hist(treats_randomizr)

data <- data.frame(values = c(treats_quaest,
                              treats_randomizr),
                   group = c(rep("survey company", 1000),
                             rep("randomizr", 1000)))

#Creating Appendix Figure A3
pdf(paste0(home, "/outputs/fig-A1.pdf"))
ggplot(data, aes(x = values, fill = group)) +
  geom_histogram(position = "identity", alpha = 0.4, bins = 50) +
  labs(x = "Number of Units assignment to Treatment", 
       y = "Count")
dev.off()

#Descriptive stats
print(xtable(rbind(c(summary(treats_quaest) %>% unname(), sd(treats_quaest)), 
       c(summary(treats_randomizr) %>% unname(), sd(treats_randomizr)))), 
      type = "latex", 
      file = paste0(home, "/outputs/tab-A4.tex"))


#Probability calculator (cited in paper)
pbinom(575, size = 1000, prob = 0.5, lower.tail = F)

#Other relevant quantities
pbinom(425, size = 1000, prob = 0.5)
pbinom(575, size = 1000, prob = 0.5, lower.tail = F) + pbinom(425, size = 1000, prob = 0.5)
dbinom(575, size=1000, prob=0.5)



